function [ output subjectinfo ] = clusterTfcPeaks( subjectinfo, cfg )
%CLUSTERTFCPEAKS

if ~exist('subjectinfo', 'var'); error('No subject information..'); end;
if ~exist('cfg','var'); cfg = clusterTfcPeaks_DefaultCfg(); end;


clusters = [];
for s = cfg.EpochSet
    
    if any(strcmp(s, fieldnames(subjectinfo.EpochData)))
        disp(['Loading epochset ''' subjectinfo.EpochData.(cell2mat(s)).epochedSetname '''...']);
        load(subjectinfo.EpochData.(cell2mat(s)).TfcDatafile);
        
        % Correctie voor tfc.times
        tfc.times = tfc.times - 5000;
        
        datatable = BuildDatatable(tfc, subjectinfo);
        X = cell2mat(datatable(2:end, 5));
        Y = cell2mat(datatable(2:end, 6));
        Z = cell2mat(datatable(2:end, 4));
        T = cell2mat(datatable(2:end, 7));
        B = datatable(2:end, 3);
        
        regions = unique(datatable(2:end, 1));
        bins = unique(Z);
        for b = bins'
            bintf = Z == b;
            if ~any(bintf); continue; end; % if no peaks in bin, skip..
            binName = B{find(bintf, 1, 'first')}; % Get the first of the list..
            for r = regions'
                region = cell2mat(r);
                regtf = strcmp(region, datatable(2:end, 1));
                tf = bintf & regtf;
                if ~any(tf); continue; end; % if no peaks in region (within bin), skip..
                dataValues = Y(tf);
                dataTimes = tfc.times(X(tf))';
                types = T(tf);
                clusters = [clusters FindClusters(dataTimes, dataValues, tfc.subject, tfc.epochs, binName, region, types)];
            end
        end
    else
        disp(['No epochset of type ''' cell2mat(s) ''' for ' subjectinfo.Name]);
    end
end

output.clusters = clusters;

end



function datatable = BuildDatatable(tfc, subjectinfo)
datatable = {'Region', 'Electrode', 'Bin', 'BinNr', 'Latency', 'Value', 'Type'};
electrodes = fieldnames(tfc.data);
for i=1:length(electrodes)
    electrode = cell2mat(electrodes(i));
    currentRow.Electrode = electrode;
    currentRow.Region = ElectrodeInfo.RegionByLabel(subjectinfo.Electrodes, electrode);
    for j=1:length(tfc.data.(electrode));
        d = tfc.data.(electrode)(j);
        % Find the bin description
        currentRow.BinNr = find(tfc.freqbins(:,1)==d.frequencyrange(1) & tfc.freqbins(:,2)==d.frequencyrange(2));
        currentRow.Bin = tfc.freqbindescriptions{currentRow.BinNr};
        for p = 1:length(d.peaklocs)
            currentRow.Latency = d.peaklocs(p);
            currentRow.Value = d.mean(d.peaklocs(p));
            currentRow.Type = d.peaktype(p);
            % Write the data to the table..
            datatable(end+1,:) = { ...
                currentRow.Region, ...
                currentRow.Electrode, ...
                currentRow.Bin, ...
                currentRow.BinNr, ...
                currentRow.Latency, ...
                currentRow.Value, ...
                currentRow.Type};
        end
    end
end
end

function [clusters] = FindClusters(times, vals, subject, epochset, bin, region, types)
% Precondition-check
if ~(size(times)==size(vals)); error('Severe problem..'); end;

% Get the cluster membership index
membership = AssignMembership(times, vals, types);
% disp(['Found ' num2str(length(unique(membership))) ' cluster(s).']);

% Build the data structure
clusters = [];
for c = unique(membership)
    cluster = [];
    members = membership==c;
    T = times(members);
    cluster.members.time = T;
    V = vals(members);
    cluster.members.value = V;
%     for i = 1:length(members)
%         cluster.members(i).time = T(i);
%         cluster.members(i).value = V(i);
%     end
    
    % Reject cluster if <2 members
    %   if length(members)<2; continue; end;
    
    % Reject cluster based on chance
    %   tcluster= max(T)-min(T);
    %   lambda=length(times)/(8000/tcluster);
    %   kcluster=length(members);
    %   chance = (exp(-lambda)*lambda^kcluster)/factorial(kcluster);
    %   if chance>0.05; continue; end;
    %     if weight<0.8; continue; end;
    %      weight = length(members)/electr; %% ratio of # of point in
    %      cluster and # of electrodes in the region
    
    % Calculate stat info
    cluster.meanTime = mean(T);
    cluster.meanValue = mean(V);
    cluster.stdTime = std(T);
    cluster.stdValue = std(V);
    
    %Calculate regression
    %     if length(V)<2;
    %      cluster.regression.coeff = NaN;
    %     else
    %      cluster.regression.coeff = regress(V-mean(V), (T-mean(T))/1000);
    %     end;
    
    % Calculate regression
    if length(V)<3; % Clusters of 1 or 2 peaks are excluded
        cluster.regression.coeff = NaN;
        cluster.regression.confidence(1)=NaN;
        cluster.regression.confidence(2)=NaN;
    else
        warning('off','stats:regress:NoConst');
        [cluster.regression.coeff, ...
            cluster.regression.confidence, ...
            R, RINT, ...
            cluster.regression.statistics] = regress(V-mean(V), (T-mean(T))/1000); %regress function always fits through origin
        warning('on','stats:regress:NoConst');
    end;
    
    %     % Calculate regression (MARIJE)
    %     if length(V)>1;
    %         fit = polyfit((T-mean(T))/1000,V-mean(V),1);
    %         cluster.regression2.polyfit = fit(1);
    %     else
    %         cluster.regression2.polyfit = NaN;
    %     end
    %     R = corrcoef(V, T);
    %     if length(V)>2
    %         C = R(1,2);
    %     else
    %         C = 0;
    %     end
    %     cluster.regression2.error = (1-abs(C))/20;
    
    % add distribution
    %    D = [times(members)' vals(members)];
    %    [n, d] = size(D);
    %   if n > d
    %       clusters(currentCluster).gmdistribution = gmdistribution.fit(D,1);
    %   else
    %       clusters(currentCluster).gmdistribution = [];
    %    end
    % add other info
    
    % Add extra info
    cluster.subject = subject;
    cluster.epochset = epochset;
    cluster.bin = bin;
    cluster.region = region;
    
    % Add the new cluster to the list
    clusters = [clusters cluster];
    
end

% minn=NaN(size(clusters));maxx=NaN(size(clusters));
% for i=1:length(clusters)
%     minn(i) = min([clusters(i).members.time]);
%     maxx(i) = max([clusters(i).members.time]);
% end
% diff = abs(maxx-minn);
% [D,I] = max(diff);
% disp(['Diff within clusters: ' num2str(D) ' in cluster ' num2str(I) '/' num2str(length(clusters))]);

end

function membership = AssignMembership(times, vals, types)

% % Simple neighbouring algorithm..
% % 1. Assign a cluster per peak
% membership = 1:length(times);
% % 2. Bubble through to find neighbours
% done = false;
% while ~done
%     OkI = true;
%     for i = 1:length(membership)
%         OkJ = true;
%         for j = 1:length(membership)
%             if i == j; continue; end; % Don't compare the same peak
%             if membership(i) == membership(j); continue; end; % Don't compare if they are already in one cluster
%             if abs(times(i) - times(j)) < 300 && sign(vals(i)-1) == sign(vals(j)-1); %TODO Should be using PeakType!!
%                 % Find all members of the second cluster, reassign all members of that cluster
%                 membership(membership==membership(j)) = membership(i);
%                 OkJ = false;
%             end
%         end
%         if ~OkJ; OkI = false; end; % If for any J memberships changed, I's will be revisited again.
%     end
%     if OkI; done = true; end;
% % Repeat 2 until no more neighbours can be joined
% end
% % 3. Reorder membership
% [B, I, membership] = unique(membership);

%% TRY BUCKET SORTING...
% For every node, select which bin should be used
% In case node fits in more than one bin, simply merge those bins
nodes = 1:length(times);
buckets = {};
for n = nodes
    % Find a bucket, or create a new one..
    candidate = false(size(buckets));
    for b = 1:length(buckets) % Loop through buckets
        bucket = buckets{b};
        for i = 1:length(bucket) % Loop through current content of bucket
             if abs(times(n) - times(bucket(i))) < 300 && types(n) == types(bucket(i)); %TODO Should be using PeakType!!
                 candidate(b) = true;
             end
        end
    end
    if any(candidate)
        % Resolve rivaling buckets
        if sum(candidate)>1
            mergeIdx = find(candidate);
            for m = 2:length(mergeIdx)
                buckets{mergeIdx(1)} = [buckets{mergeIdx(1)} buckets{mergeIdx(m)}];
                buckets{mergeIdx(m)} = [];
                candidate(mergeIdx(m)) = false;
            end
        end
        % Only one left, add to bucket
        buckets{candidate} = [buckets{candidate} n];
    else
        % Create the new bucket with node as first element
        buckets{end+1} = n;
    end
end
% Translate buckets to membership
membership = zeros(1,length(times));
for b = 1:length(buckets)
    if ~isempty(buckets{b})
        membership(buckets{b}) = max(membership) + 1;
    end
end
if any(membership==0); error('Node without membership.'); end





end

%% Old clustering algo. 
% [times Idx] = sort(times);
% vals = vals(Idx);
% membership = ones(length(times),1)*NaN;
% nextCluster = 1;
% for i=1:length(times)
%     if isnan(membership(i));
%         membership(i) = nextCluster;
%         nextCluster = nextCluster+1;
%     end
%     for j=(i+1):length(times)
%         % If the distance in time is <300 and the type of peak (Neg or Pos)
%         % is the same, they belong to the same cluster.
%         if abs(times(i) - times(j)) < 300 && sign(vals(i)-1) == sign(vals(j)-1);
%             membership(j) = membership(i);
%         end
%     end
% end

